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Abstract 

This  report  summarizes  the  research  objectives  achieved  in  this  project  during  the  period  03-01-2013 
to  02-29-2016.  Computational  techniques  have  been  developed  that  enable  large-scale  real-space 
electronic  structure  calculations  using  Kohn-Sham  density  functional  theory.  In  particular,  the 
various  components  of  the  developed  techniques  include  (i)  real-space  formulation  of  Kohn-Sham 
density- functional  theory  (DFT)  for  both  pseudopotential  and  all-electron  calculations  based  on  a 
finite-element  discretization;  (ii)  development  of  efficient  and  scalable  algorithms  for  the  solution 
of  the  Kohn-Sham  eigenvalue  problem  based  on  Chebyshev  filtering;  (iii)  development  of  reduced- 
order  scaling  techniques  by  employing  a  subspace  projection  technique  in  conjunction  with  local¬ 
ization  techniques  that  construct  a  non-orthogonal  localized  basis  spanning  the  Chebyshev  filtered 
subspace;  (iv)  leveraging  the  localized  representation  to  construct  the  quasi-continuum  reduction. 
Most  aspects  of  these  developments  have  been  numerically  implemented,  and  the  benchmark  stud¬ 
ies  have  demonstrated  that  the  developed  techniques  significantly  outperform  existing  conventional 
DFT  implementations  in  terms  of  computational  efficiency,  scaling  with  system  size,  and  parallel 
scalability  of  the  numerical  implementation.  We  believe  these  developed  techniques  and  the  numeri¬ 
cal  implementation  can  aid  as  a  platform  for  further  development  of  scalable  and  efficient  real-space 
Kohn-Sham  DFT  calculations  on  tens  of  thousands  of  atoms,  and  enable  continuing  efforts  towards 
a  seamless  bridging  of  the  quantum  and  continuum  length-scales. 


Summary  of  Research  Objectives  Achieved 


In  this  project  we  have  conducted  a  program  of  research  to  develop  computational  techniques 
that  enable  large-scale  real-space  Kohn-Sham  DFT  calculations  which  also  form  the  basis  for  the 
quasicontinuum  reduction  that  seamlessly  bridge  the  quantum  length-scale  with  the  continuum  and 
enable  an  accurate  electronic  structure  study  of  defective  crystals.  A  comprehensive  summary  of 
the  important  research  objectives  achieved  in  this  project  is  presented  below  with  references  to  the 


published  articles  (Motamarri  et  ah,  2013  Iyer  et  al.,  2014  Motamarri  &  Gavini  2014,  2016). 
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1  Real-space  formulation  of  Kohn-Sham  DFT  and  finite-element 
discretization 


1.1  Real-space  formulation 

The  ground-state  energy  of  a  materials  system  in  Kohn-Sham  DFT  is  given  by 


E(p,  R)  —  Ts(p )  +  Exc(p )  +  En(p)  +  Eext(p ,  R)  +  EZZ(R) .  (1) 


where  p  is  the  ground-state  electron-density,  R  =  {Ri,  R2,  ■  •  •  ,-Rva}  denotes  the  collection  of 
nuclear  positions  in  the  system  with  the  Ith  nucleus  having  an  atomic  number  Zj.  For  sake  of 
simplicity,  we  present  the  ideas  in  a  non-periodic  setting,  but  extension  to  periodic  and  semi-periodic 
problems  is  straightforward  (as  the  formalism  is  in  real-space)  and  has  already  been  developed  and 
implemented.  In  equation  ([Tj) ,  the  first  term,  Ts(p),  represents  the  kinetic  energy  of  non-interacting 
electrons,  the  second  term,  Exc(p),  represents  the  exchange-correlation  term  which  accounts  for 
the  quantum-mechanical  interactions  between  electrons,  and  p(x)  represents  the  electron-density 
distribution.  The  last  three  terms  of  the  energy  functional  ([!])  represent  the  classical  electrostatic 
interactions — electron-electron,  electron-nuclei  and  nuclei-nuclei  interaction  energies,  respectively, 
and  are  given  by: 


eh(p)  =  2 


P(x)p(x') 


x  —  X' 


dx  dx.'  ;  Eext(p,R) 


Zi  d(x) 
x  —  Rj 


dx  ; 


EZZ(R ) 


1  y  Z:Zj 

2  Z— ^  I  Rj  —  Jij 

i,j= 1  1  1  J 


The  non-local  (extended)  electrostatic  interactions  involve  the  kernel  1/ 1 r  —  r7|,  which  is  the  Green’s 
function  of  the  Laplace  operator.  These  electrostatic  interactions  are  thus  computed  by  taking 
recourse  to  the  solution  of  the  Poisson  equation  to  evaluate  the  total  electrostatic  potential,  </?(x,  R), 


associated  with  the  total  charge  distribution  (Motamarri  et  ah,  2012,  2013).  Finally,  the  electron- 


density  p  and  the  kinetic  energy  of  non-interacting  electrons  are  computed  from  the  self-consistent 
held  (SCF)  iteration  of  the  Kohn-Sham  eigenvalue  problem  which  results  from  the  first  variation 
of  the  DFT  energy  functional  in  equation  ([Tj) .  The  system  of  equations  corresponding  to  the 
self-consistent  formulation  are: 


V2  +  VeS(p,  R)  \  ipi  =  ei'ipi,  VeS(p,  R)  =  +  <p(p,  R)  (2a) 

N 

=  Ne\  p(x)=2^/(q)|Wx)|2,  (2b) 

i  z—1 

Na 

-—V2tp{x,R)=p(x.)-'^2zi5(x-R1)  (2c) 


7=1 


where  i/’i  denote  the  wavefunctions  or  orbitals  of  a  material  system  with  Ne  electrons  and  fi  denotes 
the  orbital  occupancy  given  by  the  finite-temperature  Fermi-Dirac  distribution.  Upon  solving  the 
above  system  of  equations  for  lowest  N  eigenfunctions,  the  kinetic  energy  is  given  by  Ts(p)  = 
Yli  pVeE-  In  the  case  of  pseudopotential  calculations,  where  only  the  wavefunctions  for  the 

valence  electrons  are  computed,  the  Coulomb  singular  Vext  is  replaced  by  a  smooth  pseudopotential 
which  provides  an  effective  potential  of  the  nucleus  and  the  core  electrons. 


1.2  Finite-element  discretization 

We  use  the  finite-element  (FE)  basis,  a  C°  continuous  piecewise  polynomial  basis,  for  discretizing 
the  Kohn-Sham  Hamiltonian.  The  FE  basis  functions  are  local  in  nature  and  result  in  discrete  sparse 
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matrices,  and,  thus,  numerical  implementations  based  on  finite-element  basis  scale  well  on  massively 
parallel  computing  platforms  (up  to  thousands  of  processors)  and  offers  a  distinct  advantage  in 
comparison  to  other  widely  used  basis  sets.  The  discretization  of  Kohn-Sham  equations  ©  using 
a  finite-element  basis  results  in  a  generalized  eigenvalue  problem  given  by 

=  e?M^i  ,  (3) 


where  H jk  denotes  the  discrete  Hamiltonian  matrix,  denotes  the  overlap  matrix  and  denotes 
jth  eigenvaiue  corresponding  to  the  eigenvector  SP  L  which  corresponds  to  the  nodal  vector  of  the 
corresponding  eigenfunction  'i/’f  (x)  on  a  given  finite-element  mesh.  The  sparse  matrices  H  is  given 

by 

ttjk  =  l  [  [viVj  (x)  .VlVfc(x)  +  fl(Ax),  fl)^(x)JVfc(x)l  dx  (4) 


where  Nj  :  1  <  j  <  M  denotes  the  finite-element  basis  functions  and  Mjk  =  f  Nj(x)Nh(x)  dx 
denotes  the  overlap  matrix.  The  key  strategies  employed  to  improve  the  computational  efficiency 
of  FE  discretization  of  Kohn-Sham  equations  are: 


A  priori  mesh  adaption:  An  a  priori  mesh  adaption  technique  is  developed  to  construct  a  close 
to  optimal  FE  discretization  of  the  Kohn-Sham  DFT  problem  by  minimizing  the  FE  approximation 
error  in  the  ground-state  energy  for  a  fixed  number  of  elements  in  the  mesh  (Motamarri  et  al. ,  2013). 
To  achieve  this,  a  formal  bound  on  the  approximation  error  is  derived,  as  a  function  of  mesh  size, 
degree  of  FE  basis  functions  and  the  electronic  fields  involved.  This  technique  enables  us  to  exploit 
the  adaptive  nature  of  FE  basis  sets  by  systematically  coarse-graining  the  finite-element  mesh  in 
order  to  reduce  the  total  number  degrees  of  freedom  and  thereby  achieve  significant  computational 
savings. 


Higher-order  spectral  finite-elements:  Higher-order  spectral  finite-elements  up  to  sixth  order 
are  employed.  Here  the  FE  basis  functions  corresponds  to  Lagrange  polynomials  interpolated 
through  an  optimal  distribution  of  nodes  (Gauss-Lobatto-Legendre  (GLL)  points),  corresponding 
to  the  roots  of  derivatives  of  Legendre  polynomials.  The  use  of  spectral  finite-elements  in  higher- 
order  FE  discretization  of  the  Kohn-Sham  eigenvalue  problem  offers  two  main  advantages.  Firstly,  a 
spectral  finite-element  basis  has  an  improved  conditioning  in  comparison  to  the  conventional  finite- 
elements  using  equispaced  nodes.  Secondly,  employing  a  GLL  quadrature  rule  in  conjunction  with 
spectral  finite-elements  results  in  a  diagonal  overlap  matrix  which  allows  for  trivial  transformation 
of  the  generalized  eigenvalue  problem  ^  to  a  standard  eigenvalue  problem  as 

H =  ejM =>  =>  H  Pt  =  (5) 

where  P,  =  H  =  M^1/2HM'1/2 


We  conducted  systematic  studies  to  assess  the  computational  efficiency  afforded  by  higher-order 
FE  discretizations  of  Kohn-Sham  DFT  (Motamarri  et  ah,  2013).  Our  benchmark  numerical  studies 
demonstrated  optimal  rates  of  convergence  for  both  pseudopotential  as  well  as  all-electron  calcula¬ 
tions  (cf.  Figures  1  &  2).  While  higher-order  discertizations  are  more  accurate  and  require  lesser 
number  of  basis  functions  to  achieve  chemical  accuracy  in  comparison  to  a  lower-order  discretiza¬ 
tion,  the  per  basis  computational  cost  is  also  higher  for  higher-order  discretizations.  Thus,  in  order 
to  assess  the  computational  efficiency  of  higher  order  finite-element  discretizations  we  conducted 
a  study  to  compare  the  CPU  times  required  to  compute  the  electronic  structure  of  benchmark 
problems  using  various  orders  of  finite-element  basis  functions.  Our  studies  indicate  that  a  4th 
order  FE  basis  function  provides  a  staggering  100-1000  fold  computational  savings  in  comparison 
to  a  linear  FE  basis  function  (cf.  Figures  3  &  4).  This  clearly  demonstrates  the  computational 
efficiency  afforded  by  higher-order  discretizations  in  electronic  structure  calculations. 
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Figure  1:  Convergence  of  error  in  ground-state  en¬ 
ergy  of  a  35-atom  Barium  cluster  for  various  orders 


Figure  2:  Convergence  of  error  in  ground-state 
energy  of  a  methane  molecule  for  various  orders 


of  finite-elements  (pseudopotential)  (Motamarri  et  of  finite-elements  (all-electron)  (Motamarri  et  al. 


al. , 

2013 

)• 

2013) 

1.3  Solution  of  the  Kohn-Sham  eigenvalue  problem  via  Chebyshev  filtering: 

The  use  of  spectral  elements  allows  us  to  trivially  transform  the  generalized  eigenvalue  problem 
to  a  standard  eigenvalue  problem.  Noting  that  the  key  quantity  of  interest  in  density  functional 
theory — electron-density — can  be  computed  from  the  eigenspace  of  the  occupied  states,  we  em¬ 
ploy  the  Chebyshev  filtering  technique  (ChFSI-FE)  to  compute  an  approximation  to  the  occupied 
eigenspace.  We  refer  to  Motamarri  et  al.  (2013);  Motamarri  &  Gavini  (2014)  for  the  specific 


details  of  our  implementation.  In  brief,  the  ChFSI  technique  represents  a  modified  power  itera¬ 
tion  in  every  SCF  iteration  that  exploits  the  rapid  growth  of  Chebyshev  polynomials  in  (— oo,  —1) 


to  magnify  the  relevant  eigenspectrum  (Zhou  et  al.,  2006),  thus  providing  a  good  subspace  for 


the  solution  of  the  Kohn-Sham  eigenvalue  problem.  Upon  computing  the  Chebyshev  filtered  sub¬ 
space,  we  orthogonalize  the  Chebyshev  filtered  vectors  and  project  the  discrete  Hamiltonian  onto 
the  subspace  and  compute  the  eigenvalue  and  eigenvectors  of  the  subspace  projected  Hamiltonian. 
We  evaluated  the  efficiency  of  the  ChFSI-FE  approach  with  conventional  approaches  of  solving 
the  generalized  Kohn-Sham  eigenvalue  problem  using  Jacobi-Davidson  method  and  the  standard 
eigenvalue  problem  using  the  Krylov  Schur  method.  Table  1  shows  this  comparison  demonstrat¬ 
ing  the  significant  computational  efficiency  afforded  by  ChFSI-FE  for  both  pseudopotential  and 
all-electron  calculations. 

Table  1:  Comparison  of  generalized  eigenvalue  problem,  standard  eigenvalue  problem,  and  the  Chebyshev 
filtering  technique.  Benchmark  problems  include  an  all-electron  calculation  of  16  atom  graphene  sheet  and 
a  pseudopotential  calculation  on  a  172  atom  aluminum  cluster  (Motamarri  et  al. , |2013). 


Element  Type 

DOFs 

Problem  Type 

N 

Time  (GHEP) 

Time  (SHEP) 

ChFSI-FE 

HEX125SPECT 

1,368,801 

graphene 

96 

1786  CPU-hrs 

150  CPU-hrs 

12.5 

HEX343SPECT 

2,808,385 

Al  3  x  3  x  3  cluster 

516 

2084  CPU-hrs 

80  CPU-hrs 

13 

We  also  conducted  benchmark  studies  comparing  the  developed  approach  using  ChFSI-FE  with 
ABINIT  software,  which  is  a  well  developed  plane-wave  based  method  for  solving  the  ground-state 
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Figure  3:  Relative  performance  of  various  orders 
of  finite-elements  measured  in  terms  of  computa¬ 
tional  time  for  pseudopotential  calculation  of  35- 
atom  Barium  cluster  (Motamarri  et  al.  2013). 


Figure  4:  Relative  performance  of  various  orders  of 
finite-elements  measured  in  terms  of  computational 
time  for  the  all-electron  calculation  of  methane 
molecule  (Motamarri  et  al.  2013). 


Kohn-Sham  DFT  problem.  Table  2  shows  these  results  which  demonstrates  a  significant  outperfor- 
mance  by  the  developed  approach.  Using  the  developed  real-space  approach,  we  have  demonstrated 
pseudopotential  calculations  on  systems  as  large  as  a  few  thousand  atoms  and  all-electron  calcula¬ 
tions  on  systems  as  large  as  a  few  thousand  electrons  using  the  modest  computational  resources  at 
our  disposal.  Figures  5  &;  6  show  the  electron  density  contours  from  all-electron  calculations  of  a 
graphene  sheet  contating  100  atoms  and  Tris  (bipyridine)  ruthenium  complex,  an  organo-metallic 
complex.  Further,  given  the  good  scaling  exhibited  by  the  developed  implementation  on  parallel 
computing  architectures  (cf.  Figure  7),  systems  beyond  the  reach  of  conventional  plane- wave  based 
implementations  are  now  accessible. 

Table  2:  Comparison  of  the  developed  real-space  DFT  code  with  ABINIT  for  an  Aluminum  nanocluster 
containing  666  atoms  (pseudopotential  calculation  with  Troullier  Martins  psedopotential) . 


Type  of  basis  set 

Energy  (eV) 
per  atom 

Abs.  error  (eV) 
per  atom 

Rel.  error 

Time  (CPU-hrs) 

ABINIT  (cut-off  20  Ha-, 
1,093,421  plane  waves) 

-56.8739 

0.0043 

7.5  x  10"5 

8640 

Real-space  (HEX343SPECT 
FE  basis;  4,363,621  nodes) 

-56.8765 

0.0017 

3  x  10-5 

1132 

Reduced  order  scaling  using  localization  ideas  and  quasi-continuum 
reduction: 


The  real-space  ChFSI-FE  formulation  discussed  in  the  previous  section,  which  computes  the  canoni¬ 
cal  eigenfunctions  and  eigenvalues,  asymptotically  scales  cubically  with  number  of  atoms  (due  to  the 
cubic  scaling  of  diagonalizing  a  matrix).  This  becomes  prohibitively  expensive  for  large-scale  DFT 
calculations,  and  hence,  in  order  to  reduce  the  computational  complexity,  we  developed  a  subspace 
projection  technique  (Motamarri  &  Gavini,  2014)  that  has  reduced  order  scaling.  The  formulation 
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Figure  5:  Electron-density  contours  from  an  all-electron 


DFT  calculation  of  a  100  atom  graphene  sheet  (Motamarri 


et  al.  2013). 


Figure  6:  Isocontours  of  the  electron 
density  from  an  all-electron  calcula¬ 
tion  of  Tris  (bipyridine)  ruthenium 
complex  (Motamarri  et  al.  2013). 


Figure  7:  Strong  scaling  study  demonstrating  over  90%  efficiency  at  768  processors.  Benchmark 
system  is  3  x  3  x  3  aluminum  cluster  (172  atom)  with  a  FE  discretization  containing  4  million 
basis  functions. 
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provides  a  single  unified  framework  to  handle  both  insulating  and  metallic  systems.  Further,  both 
pseudopotential  as  well  as  all-electron  calculations  can  be  conducted  using  the  proposed  method¬ 
ology.  The  key  ideas  involved  in  the  method  include:  (i)  expressing  the  Kohn-Sham  Hamiltonian 
and  the  corresponding  wave-functions  in  Lowdin  orthonormalized  finite-element  basis  constructed 
using  spectral  finite-elements  in  conjunction  with  GLL  quadrature  rules;  (ii)  employing  a  Cheby- 
shev  filtering  approach  in  each  SCF  iteration  for  an  efficient  computation  of  the  eigenspace;  (iii) 


use  a  localization  procedure  (Garcia-Cervera  et  ah, 

2009) 

to  compute  atom-centered  nonorthogonal 

localized  basis  (localized  wavefunctions)  spanning  the  Chebyshev  filtered  subspace;  (iv)  comput¬ 
ing  the  projection  of  the  Hamiltonian  into  the  non-orthogonal  localized  basis  and  using  a  Fermi- 

operator  (Baer  &  Head-Gordon,  1997)  expansion 

of  the  projected  Hamiltonian  to  compute  the 
nsity  and  the  band  energy.  The  computational 
refunctions  with  a  compact  support  is  shown  to 

relevant  quantities  like  density  matrix,  electron  de 
complexity  of  this  algorithm  for  well-localized  wav 

be  linear  -  O(N).  We  refer  to  Motamarri  &  Gavini 

(2014 

)  for  details  of  this  algorithm,  and  present 

here  some  benchmark  results  that  demonstrates  the  significant  computational  efficiency  afforded 
by  this  approach.  Figure  8  shows  the  computational  time  per  SCF  for  the  developed  subspace  pro¬ 
jection  method  in  comparison  with  ChFSI-FE  for  benchmark  problems  involving  a  series  aluminum 
nano-clusters  and  alkane  chains.  As  is  evident  the  developed  localization  approach  provides  close 
to  ten-fold  computational  savings  and  enables  consideration  of  systems  approaching  close  to  10,000 
atoms.  Further,  we  note  that  the  observed  scaling  for  the  alkane  chains  (system  with  band  gap) 
is  close  to  linear  scaling  with  system  size,  and  the  scaling  for  aluminum  clusters  (representative 
of  a  metallic  system)  is  ~  0(N1'5),  which  is  a  significant  improvement  over  the  cubic-scaling  for 
conventional  solution  methods. 


(b)  Insulating  system 


Figure  8:  Average  computational  time  (CPU-Hrs)  per  SCF  iteration  for  the  subspace  projection  method 
and  ChFSI-FE  for  a  pseudopotential  calculation.  Case  study:  (a)  Aluminum  nano-clusters(upto  3430  atoms) 


(b)  Alkane  chains  (up  to  7000  atoms)  (Motamarri  &  Gavini,  2014 1 . 


While  providing  reduced  order  scaling,  the  developed  localization  technique  is  also  crucial  to 
the  quasi- continuum  reduction  of  Kohn-Sham  DFT.  Using  this  localized  representation,  the  wave- 
functions  are  decomposed  into  a  predictor  and  corrector,  where  the  predictor  is  computed  using  the 
Cauchy  Born  hypothesis  (along  the  similar  lines  as  the  quasi-continuum  reduction  for  orbital-free 
DFT  (Gavini  et  al.  2007 ;  [Radhakrishnan  fc  Gavini  2010)).  In  regions  of  smooth  deformations,  such 
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as  away  from  the  defect-core,  the  predictor  provides  an  effective  representation  of  the  wavefunc- 
tion,  whereas  in  regions  close  to  the  defect-core  the  corrector  provides  corrections  to  the  predictor. 
Given  the  nature  of  the  corrector  fields,  they  can  be  efficiently  represented  by  using  an  adaptive 
discretization,  with  more  resolution  near  the  defect-core  and  coarse-graining  away  from  it.  We  are 
currently  implementing  the  quasi-continuum  reduction  for  Kohn-Sham  DFT  using  the  localization 
technique  in  our  framework. 


3  Application  studies: 


We  investigated  the  energetics  of  point  defects — monovacancy,  self-interstitials  (tetrahedral,  octa¬ 
hedral  and  dumbbell)  and  Frenkel  pairs — in  aluminum  under  tensile  hydrostatic  stresses  in  order  to 
understanding  the  mechanisms  leading  to  spall  failure  in  metals.  We  find  that  that  the  defect  core¬ 
energy  of  a  vacancy,  which  is  governed  solely  by  the  electronic-structure  at  the  core,  significantly 
depends  on  the  macroscopic  hydrostatic  stress,  and  that  this  constitutes  an  important  contribution 
to  the  formation  enthalpy,  especially  in  the  regime  of  extreme  tensile  hydrostatic  stresses.  This 
finding  is  in  contrast  to  widely  used  elastic  formulations  of  point  defects  based  on  formation  vol¬ 
ume  that  ignore  the  defect  core-energy  contribution.  The  formation  enthalpies  of  all  point  defects 
considered  in  this  study  monotonically  decrease  with  increasing  tensile  hydrostatic  stress  (cf.  Fig¬ 
ure  9).  Further,  we  find  that  the  formation  enthalpies  of  vacancies  and  Frenkel  pairs  are  negative 
beyond  critical  tensile  hydrostatic  stresses  (9  GPa  for  vacancies  and  12  GPa  for  Frenkel  pairs), 
which  suggests  a  spontaneous  nucleation  of  these  point  defects  and  this  has  important  implications 
to  nucleation  mechanisms  leading  to  spall  failure.  In  particular,  these  findings  suggest  two  possible 
defect  nucleation  mechanisms  leading  to  spall  failure:  (i)  a  heterogeneous  nucleation  of  vacancies 
from  defect  sources;  (ii)  a  homogeneous  nucleation  of  Frenkel  pairs  at  higher  hydrostatic  stresses. 
We  refer  to  (Iyer  et  al.,  2014)  for  further  details  on  the  calculations  and  the  findings. 
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-  Self-interstitial  Formation  Enthalpy  (eV) 


(a)  Vacancy 


(b)  Self  interstitial 


(c)  Frenkel  pair 


Figure  9:  Formation  enthalpies  of  point  defects  as  a  function  of  tensile  hydrostatic  stress  (Iyer  et  al.  2014 1 . 
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